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Abstract. The Cauchy problem for the Korteweg de Vries (KdV) equation with small 
dispersion of order e^, e <C 1, is characterized by the appearance of a zone of rapid mod- 
l/^ [ ulated oscillations of wave-length of order e. These oscillations are approximately de- 

■ scribed by the elliptic solution of KdV where the amplitude, wave-number and frequency 

I are not constant but evolve according to the Whitham equations. In this manuscript 

we give a quantitative analysis of the discrepancy between the numerical solution of the 
^ I KdV equation in the small dispersion limit and the corresponding approximate solution 

^ O ■ for values of e between 10~^ and 10~^. The numerical results are compatible with a 

difference of order e within the 'interior' of the Whitham oscillatory zone, of order 
• at the left boundary outside the Whitham zone and of order ^/e at the right boundary 

outside the Whitham zone. 
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1. Introduction 

The purpose of this manuscript is the quantitative numerical comparison of the solution 



O ■ of the Cauchy problem for the Korteweg de Vries equation (KdV) 

Q^. (1-1) ut + 6uu^ + e^U:^^^ = 0, u{x,0) = Uo{x), 

I in the small dispersion limit (small e) and the asymptotic formula obtained in the works of 

^ ' Lax and Levermore |22j, Venakides |E2j and Deift, Venakides and Zhou \^ which describes 

Vh ■ the solution of the above Cauchy problem at the leading order as e — 0. We study initial 

>■ ! data with a negative hump and with a single minimum value at x = normalized to 

I —1. The solution of the Cauchy problem for the KdV equation is characterized by the 

^ • appearance of a zone of fast oscillations of wave-length of order e, see e.g. Fig. ^ These 

oscillations were called by Gurevich and Pitaevski dispersive shock waves fl^ . 
Following the work of [23; IHZ] and [S], the rigorous theoretical description of the small 
dispersion limit of the KdV equation is the following. Let us define 

(1.2) u{x,t) = \imu{x,t,e). 

1) for < t < tc, where is a critical time, the solution u{x,t,e) of the KdV Cauchy 
problem is approximated, for small e, by the limit u{x,t) which solves the Hopf equation 

(1.3) Ut + 6uUx = 0. 
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Figure 1. The numerical solution of the KdV equation at different times 
for the initial data uo{x) = — 1/ cosh^x and e = 10~^'^. 

Here tc is the time when the first point of gradient catastrophe appears in the solution 
(1.4) u{x,t) = uoiO, x = 6tuo{0 + ^, 

of the Hopf equation. From the above, the time tc of gradient catastrophe can be evaluated 
from the relation 

1 

tr = min 



2) After the time of gradient catastrophe, the solution of the KdV equation is characterized 
by the appearance of an interval of rapid modulated oscillations. According to the Lax- 
Levermore theory, the interval [x~ (t) , (t)] of the oscillatory zone is independent of 
e. Here x~{t) and x~^{t) are determined from the initial data and satisfy the condition 
x~(tc) = x^{tc) = Xc where Xc is the x-coordinate of the point of gradient catastrophe 
of the Hopf solution. Outside the interval [x~ (t) , x~^ (t)] the leading order asymptotics 
of u{x,t,e) as e ^ is described by the solution of the Hopf equation ()1.4|) . Inside the 
interval [x~ (t) , x'^ (t)] the solution u{x, t, e) is approximately described, for small e, by the 
elliptic solution of KdV [HI, j23,EII,|HI 

.2 ,„„./v/A^, 



1.5) uix, t,e)^u + 2e'— log 6 ^^^^^^ " ^Pi + f32 + Ps) - q]; T 
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where now u = u{x, t) takes the form 

(1.6) u = (3i+ (32 + (3z + 2a, 



(1-7) a = -(3i + {(3i- (J'i)——, T = t——, s - 



K{sy K{s)' (3, -(3s 

with K{s) and E{s) the complete elhptic integrals of the first and second kind, K'{s) = 
K{y/l — s^); 6 is the Jacobi elliptic theta function defined by the Fourier series 

neZ 

For constant values of the (3i the formula is an exact solution of KdV well known in 
the theory of finite gap integration |2Qj, [U^- On the contrary, in the description of the 
leading order asymptotics of u{x, t, e) as e — >• 0, the quantities Pi depend on x and t and 
evolve according to the Whitham equations |SHI 

■TrA + Vij^Pi = 0, ^ = 1,2,3, 
ot ox 

where the speeds Vi are given by the formula 

(1.8) .. = 4lM^Z^ + 2(A+/5. + /?3), 

with a as in p.7|) . Lax and Levermore first derived, in the oscillatory zone, the expression 
()l.(i|l for u = u{x,t) which clearly does not satisfy the Hopf equation. The theta function 
formula ()1.5|1 for the leading order asymptotics of u{x,t,e) as e — 0, was derived in 
the work of Venakides and the phase q = q{Pi, P2, Ps) was derived in the work of Deift, 
Venakides and Zhou [8^ for the case of pure radiation initial data, namely initial data 
with no solitons. However, we verify numerically that their formula holds also for initial 
data with point spectrum. We give a formula for q which looks different but which is 
equivalent to the one in 



where f-{y) is the inverse function of the decreasing part of the initial data. The above 
formula holds till some time T > tc. For later times the formula has to be slightly modified 
(see formula ()2.14|) ). The function q = q{Pi, P2, Ps) is symmetric with respect to Pi, P2 
and Ps, and satisfies a linear over-determined system of Euler-Poisson-Darboux type. It 
has been introduced in the work of Fei-Ran Tian jH^ to study the Cauchy problem for 
the Whitham equations after the first breaking of the Hopf solution. A formula for the 
phase in the multi-bump case is derived in [T3] . 

3) Fei-Ran Tian proved that the description in 1) and 2) is generic for some time after 
the time tc of gradient catastrophe |S2]. 

In this paper we compare numerically, outside and inside the oscillatory region, the as- 
ymptotic formulas ()1.4j) and ()1.5j) with the solution of the KdV equation, see Fig. |21for a 
plot of the two solutions for e = 0.1. Our main observations are the following. 
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Figure 2. The blue line shows the solution to the KdV equation for e = 0.1 
for the initial data uo{x) = —1/ cosh^ x at the time t = 0.4, the green line 
the asymptotic solutions p.4|l and 



(1) The oscillations of the numerical solution of the KdV equation start at a time 
t < tc- Indeed for the initial data Uq{x) = — 1/ cosh^ x the critical time of the Hopf 
solution is given by tc = a/3/8 ~ 0.2165 and it can bee seen from Fig. IHlthat the 
KdV solution has already developed two oscillations at this time. 

(2) The oscillatory interval of the KdV solution is bigger than the oscillatory interval 
[x~ (t) , x^ (t)] described by the leading order asymptotics given by formula p.Sj) . 
From the numerical simulation it can be seen that the KdV oscillatory interval is 
shrinking, as e ^ 0, to the oscillatory interval [x~{t), x~^{t)] (see Fig. [Hand Fig.|3J. 
Let us define '^topf '■— -^i^topf/^^ ~ 1) where x^^^j are the values of x for which 
the absolute value of the difference between the KdV solution and the asymptotic 
solution ()1.4j] are smaller than some fixed value (we take it to be 10~^ here) for all 
^x > x^opf Then we numerically obtain that ^^opf ^ with standard error 
0.028 for the exponent. However, this result has to be taken with care in view of 
the low number of points and the arbitrariness in the definition of the zone. The 
zone ^hopf clearly shrinking with e but the dependence on e does not seem to 
be described by a power law. 

(3) At the boundaries of the oscillatory region, the leading order asymptotics described 
by formula ()1.5|1 matches C°, but not C^, the solution of the Hopf equation ()1.4|1 
(see Fig. ini). We prove this statement analytically in Theorem 12.21 

(4) The difference between the KdV solution and the asymptotic solution is decreasing 
with e. To define an error, we take the maximum of the absolute value of the differ- 
ence between the solutions close to the center of the Whitham zone [x~{t),x~^{t)]. 
We find that this error is of order e (more precisely e^'^^^ with standard error 
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Figure 3. The blue line is the solution of the KdV equation for the initial 
data Uo{x) = —1/ cosh^x and e = 10~^, and the purple line is the corre- 
sponding leading order asymptotics given by formulas ()1.4|) and ()1.5|) . The 
plots are given for different times near the point of gradient catastrophe 
{xc,tc) of the Hopf solution. Here Xc ^ —1.524, tc ^ 0.216. 



0.05). At the left boundary of the oscillatory zone, and for x < x (t) we obtain 
that the error is decreasing like e°'^^ with standard error 0.025 which is roughly es. 
The error for x > x~^{t) is decreasing like e'^^^ with standard error 0.017, which is 
roughly y/e. 

The manuscript is organized as follows: in the second section we give a brief analytical 
description of the Lax-Levermore theory and the Whitham equations. In the third and 
fourth section we give a description of the numerical algorithm used to solve the KdV 
equations and the Whitham equations, respectively. Readers not interested in the nu- 
merical details can skip this part. In the fifth section we present a detailed comparison 
of the numerical KdV solution in the limit of small dispersion and the corresponding as- 
ymptotic solution. In the last section we add some concluding remarks. We concentrate 
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Figure 4. The numerical solution of the KdV equation for the initial data 
uo{x) = — 1/ cosh^x, for different values of e and for fixed time t = 0.4. 
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(a) KdV solution 



(b) Asymptotic solution 



Figure 5. In (a) the numerical solution of KdV is plotted for t = 0.4 and 
e = 10"^. In (b) the asymptotic formula p.5|l and p.4|l is plotted for the 
same values of t and e. 
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Figure 6. The difference between the plots (a) and (b) of Fig. HJ the 
Whitham zone is shown in blue, the exterior of this zone in green. 

on the explicit example of initial data Uo{x) = —1/ cosh^x. The codes we developed are, 
however, suitable for any initial data with a single hump, which are rapidly decreasing at 
infinity. 

2. Lax Levermore theory and Whitham equations 

In this section we sketch the Lax-Levermore theory and its connection with the Whitham 
equations. We consider initial data with one negative hump that tends to zero fast enough 
as X =— s> ±oo. It is well known that the solution of the KdV equation can be obtained 
by the inverse scattering method. Through this explicit expression of the KdV solution. 
Lax and Levermore [29j for positive hump initial data and Venakides ^^6^ for negative 
hump initial data, managed to find the limit of the solution of as e — >■ via a vari- 
ational problem. The KdV zero-dispersion limit exists in the distribution sense and can 
be determined as follows: 

(2.1) u{x,t) = d-\imu{x,t,e) = 2dlQ{x,t) - 1, 
with 

(2.2) Qix,t) = ini Q{ij,x,t) 

where S is the set of all Lebesgue measurable function with support in [minuo(x), 0]. For 
our numerical purposes we do not use the Lax-Levermore- Venakides functional Qlip, x, t) 
and for this reason we omit the explicit formula but recall only its properties. For each 
fixed {x,t) the functional Q{ip,x,t) assumes its minimum at exactly one element of the 
set S denoted by '^*(a:, t), moreover u{x,t) can be expressed in terms of the endpoints of 
the support of ip*{x,t). 

To describe the support of the minimizer ip*{x,t), Lax and Levermore consider the 
motion of the initial curve where each point of the curve has a different speed. At time 



E=10"^ 
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t = the curve is u = Uo{x) and the support of the minimizer is [uo{x), 0]. At later times 
the curve will, in general, be given by a multivalued function in the {x, u) plane with an 
odd number of branches. For < t < tc, where is the time of gradient catastrophe for 
the Hopf equation, the curve is given by the solution u{x,t) of the Hopf equation ()1.3p 
and the support of the lax-Levermore-Venakides minimizer is [u{x, t), 0]. The limit u{x, t) 
takes the form 

u{x, t) = u{x, t). 

Soon after the time of gradient catastrophe, the evolving curve is given by the multivalued 
function parameterized by three branches > > P2{x,t) > /?3(x, t) > —1, and the 

support of the minimizer is [Ps, P2] U [Pi, 0] (see Fig. IT3|) . The limit u{x,t) is expressed by 
the formula 

■u(x, t) = (3i{x, t) + [32{x, t) + P^{x, t) + t), 

where a is defined in fll.7|) . In this case u{x,t) does not satisfy the Hopf equation ()1.3p 
but a new set of equations enter the game. Indeed the Pt as functions of x and t satisfy 
the Whitham equations 



(2.3) + /5, = 0, ^ = 1,2,3, 
where 

(2.4) v. = A +2)^/?.. 

^* fc=i 

The equations ()2.3p are also called one-phase Whitham equations to distinguish them from 
the n-phase Whitham equations |T3j which are derived when the support of the minimizer 
consists of 2n + 1 intervals. In terms of the quantities Pi,P2 and P^ the approximate 
solution of u{x, t, e) as e — > is given by the formula jSZj, [H| 



u{x, t, e) ~ + /?2 + + 2« + 2e^-^ hgO ( ^f'/' [x - 2t(A + P2 + Ps) - (j)]; T 



dx^ \ 2eK{s) 

with T and K{s) as in p.7|) . The phase = (p{Pi, P2, P3) is determined by the Deift- 
Venakides-Zhou formula [8J as follows. Let us introduce the functions p(A) and r(A) which 
are the semiclassical approximation of the reflection and transmission coefficients for the 
Schrodinger operator — e^(9^^ — Uo{x), 

x+{X) 



(2.5) p(A) = v^x+(A)+ / [V^-^/uo{y)-X\dy, r(A) = / ^/X-uoiy)dy, 

Jx+{X) Jx-{\) 

where x±{\) are the solutions of Uq{x±{X)) = A. Then the phase takes the form 

1 /• p(A)c^A 1 r tTi\)d\ 

^ Jhui, v/(A-A)(A-/32)(A-/33) ^ y[o,,]\u/, V(A-/3i)(A-/32)(A-/53)' 

where Ji = [Pi, 0], I2 = [Ps, P2] and r] = —1 when t < T while rj = [3^ when t > T where 
T is the time when (3^ first reaches the minimum value P^ = —1. 

The problem of determining the small dispersion limit of KdV is reduced to solv- 
ing the Lax-Levermore-Venakides variational problem ()2.2j) or the Whitham equations 
()2.3|) . In jSI] McLaughlin and Strain obtain numerically the weak limit u by solving the 
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Lax-Levermore-Venakides variational problem. In this manuscript we follow the latter 
approach, and we solve numerically the Whitham equations ()2.3p . thus determining the 
support of the Lax-Levermore minimizer. The numerical solution of the Whitham equa- 
tions has been already implemented in the early works of Gurevich and Pitaevski jTHI and 
Avilov and Novikov [3^ for step-like and cubic initial data. In this manuscript we imple- 
ment a code which is suitable to rapidly decreasing initial data with an arbitrary single 
hump. As a concrete example we study the case of the initial data uo{x) = — 1/ cosh^x 
in detail. 

2.1. The Cauchy problem for the Whitham equations. In this subsection we 
mainly follow the work of Fei-Ran Tian j^. The Whitham equations are a sys- 
tem of hyperbolic PDEs defined for Pi > P2 > Ps- Using the properties of the elliptic 
integrals 

(2.7) K(s) = ^ (1 + £ + + 0(s')) , E(s) = I (1 - £ - 1.= + 0(.= 



and 

(2.8) E(s)^ 1 + 1(1 



• 16 
log 



1 16 



1 - s 

we find that the speeds Vi defined in ()2.4j) have the following behavior 

1) at p2 = Pi 

viiPuPi, P3) = V2{Pi, Pi,P3) = 4/3i + 2/53 
v3{Pi,Pi,Ps) = QP3; 

2) at p2 = p3 

Vi{PuP3,P3)=QPl 

V2{Pl,P3,P3) = V3{Pi,P3,P3) = 12/3i - 6/^3; 

The initial value problem for the Whitham equations is to determine the solution of ()2.3j) 
from the following boundary conditions: 

a) Leading edge: 

Pi = the Hopf solution ()1.3p 

^^■^^ P2 = P3, 

b) Trailing edge: 

P2 = Pi 

^ ■ ^ p3= the Hopf solution (|r3|l . 

The Whitham equations are a set of quasi-linear hyperbolic PDEs jSH] that can be inte- 
grated by a generalization of the method of characteristics. Dubrovin and Novikov [T2] 
developed a geometric-Hamiltonian theory of the Whitham equations ()2.3|1 . Using this 
theory, Tsarev [SS] was able to integrate the equations through the so called hodograph 
transform, which generalizes the method of characteristics, and which gives the solution 
of ()2.3|) in the implicit form 

(2.11) X = Vit + Wi, i = 1,2,3, 
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where the Vi are defined in ()2.4p and the Wi = /?2, Ps) are obtained from an algebro- 

geometric procedure |^ by the formula ^2 

(2.12) '^^ = l(^^^-^J2f^'^^.+1^ z = 1,2,3. 

The function q is defined by 

(,13) A) r *<^.^-<^'-^- " " 



,2 



In the above formula f-{y) is the inverse function of the decreasing part of the initial 
data Uo{x). The above formula for q{Pi, P2, P3) is valid as long as Pi > P2 > P3 > — 1- 
When P3 reaches the minimum value —1 and passes over the negative hump, then it is 
necessary to take into account also the increasing part of the initial data /+ in formula 
()2.13|) . We denote by T this time. For t > T > tc we introduce the variable X3 defined 
by uoiX^) = Ps which is still monotonous. For values of X3 beyond the hump, namely 
X3 > 0, we have to substitute ()2.13|1 by the formula |^ 
(2.14) 



2vr7fe V(A-A)(A-/32)(A-/33) Vife J-iV^. 

Clearly the Wi are constructed in such a way that the matching conditions ()2.10|) and 
(I2.9|l are satisfied. Indeed 

at the trailing edge: Pi = P2 

f^Ps), fort<r 



WliPuPuPs) = W2iPl,Pl,p3), W3{pi,pi,p3) 



U{P3), f0Tt>T 



and at the leading edge: P2 = P3 

W2{Pl,P3,P3)=W3{Pl,P3,P3), Wi{Pi, P3, P3) = f-{Pl) ■ 

At the leading and trailing edge the system 1)2.111) becomes degenerate. To avoid degen- 
eracy we rewrite the system ()2.1H) in the equivalent form 

1 



(2.15) 



iPi - P2)K{s 

V3t + W3 = X 
1 



-[{Vi - V2)t + Wi~W2] =0 



[(^^2 - V3)t + W2- W3] = 0. 



In the limit P2 = P3 the system ()2.15|1 reduces to the system jHSj, [H] 

(QPit + f4Pi)-x = 

(2.16) I $(/?3,/?i) + 6t = 

[dpMP3,Pi) = 

where 
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The formula ()2.16|) has been obtained by using for the system ()2.15p the expansion ()2.7|) of 
the elhptic functions and the properties of the function q = q{Pi, P2, Ps) which is symmetric 
with respect to Pi, P2 and P3. The solution of ()2.16p yields x, Pi and P3 as a function of t. 
In particular x = x^{t) defines the left boundary of the oscillatory zone. In order to get 
an approximation near Xc of the function x = x~{t) we make a Taylor expansion of the 
system ()2.1(jj) near P3 = Pi = Uc, x = Xc and t = tc- For this purpose we first solve the 
last equation of ()2.1(j|l for P3 = Ps^Pi) and enter the second equation of ()2.1fjj) with this 
solution. Using the Taylor expansion near the point of gradient catastrophe {xc,tc,Uc), 
we obtain 

% - ~ 6{Pi - u,){t - Q + 6uc{t - Q + lfl'{u,){Pi - 





^HP.P.: 



OPi 



dpsdpi dpi dp( 



Using the identities 



/3i=/33=«c 



and 



we obtain 



92 



dPl 



m?..Pi) 



dp^P_ 



I3i=f33=uc 



MPs. Pi 



I3i=f33=uc 



/3i=/33=«c 



—f"'iu, 

^^J^K c 



a; - X, ^ 6{Pi - u,){t - Q + Qu,{t - Q + - u,f 

o 

Substituting the value of Pi — Uc from the second equation in the first, we finally obtain 

36^2 



(2.18) 



Xappit) -Xc + QUc{t - tc) 



--{t-Q 



which is the semi-cubic law obtained in JHl for cubic initial data. 
Using (Q, the limit P2 = Pi of the system ^TW} is [32], [HI 



:2.19) 



mPi,P3) + 6t = 

/3i 



/ ' ^/X-P3mX,P3)+Qt]dX = 



[ -x + 6t/33 + /-(/33) = 0, 



where the function $(A, /^s) has been defined in ()2.17j) . The solution of ()2.19j) as a function 
of t defines x = x^{t), which gives the right boundary of the oscillatory zone. As for the 
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leading edge, we can obtain the approximate behavior of the function x^{t) near the point 
of gradient catastrophe (see Fig. [7j) 



4vl0 3 

[2.20) xtppit) ~ + 6uc{t - Q + —j==={t - Q-K 




Figure 7. The growth of the Whitham zone in the {x,t) plane for the 
initial data —1/ cosh^x (blue) and for cubic initial data ()2.18p and ()2.20|) 
(green) . 



In the system (fU^ . the formula ipTTj) for $(/3i,/53) holds till the time T > when 
(3^ = —1 which can lead to /^a being non- monotonous. This time is given by the solution 
of the equations 



(2.21) 



( X + 6T = 



VA + 1[<I>(A, -1) + QT]dX = 0. 



-1 



For t > T the function <I>(A, Ps) appearing in the system ()2.19|) has to be modified to the 
form 



f2.22l 



$(A,/33 



2VX^3 




dyf'M , f'dyf'_iy) 



v/A^ 



+ 



v/A^ 



The boundary conditions ()2.9j) and ()2.1()j) guarantee that the solution of the Whitham 
equations is attached in a way in the (x, u) plane to the solution of the Hopf equation 
(see FiglHl). Namely when P2 = Ps, the derivative dxPi is continuously attached to dxu{x, t) 
where u{x,t) is the solution of the Hopf equation. The same holds in the case /?i = ^2- 
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The partial derivatives dx(3i, i = 1, 2, 3 have been obtained in 

a + A 



and take the form 



(2.23) 



dxl3i 



where Q = Q{l3i, [32, /S^) is given by 



and dxPs in the hmit 
(5, we 



i=i 

To simphfy the calculation we study d^Pi in the limit P2 
P2 Pi before passes over the negative hump. Fixing P2 = v + 6, = v 
obtain for S ^ 0, with the expansion ()2.7|) . the relation 

(2.24) 4ft = _i_ + OW. 

Relation ()2.24p shows that, in the limit 5 — > 0, the derivative d^Pi is converging to dxu{x, t) 
where u{x, t) solves the Hopf equation. Fixing [32 = v — 5,[3i = v + 5 and evaluating dxP^ 
in the limit 5 — 0, one obtains in the same way at the trailing edge 

(2.25) a,A = ^^^^ + 0(l/log^), 

which shows that (9^/53 — > d^u where u solves the Hopf equation. The solvability of the 




Figure 8. The blue line is the solution of the Whitham equations j3i > 
P2 > Ps, the green line the solution u{x,t) of the Hopf equation ()1.3j) for 
the initial data ^0(2;) = —1/ cosh^x and the time t = 0.4. The Whitham 
solution is attached to the Hopf solution. 



equations ()2.1H) and the systems ()2.1fi|l and ()2.19|1 is guaranteed by the following theorem 



14 T. GRAVA AND C. KLEIN 

Theorem 2.1. Consider smooth initial data Uq{x) with a single negative hump and sup- 
pose that Uq{x) reaches its only minimum at x = where Uq{0) = — 1. // the inverse 
function f-{u) of the decreasing part of Uq{x) satisfies the conditions 

fliu*) = 0, r(w)<0, u^u* 

where u* is the inflection point of f^{u), then the Whitham equation \2.'Jji have a unique 
solution (3i{x,t) > [32{x,t) > (3^{x,t) for t > tc till a short time after the time T. 

In the foUowing we show that the phase (f) defined in fl2.(ij) can be expressed in terms 
of q{Pi, P2, P3) defined in ()2.1Hj) or ()2.14j) . Since this quantity is aheady numerically 
implemented in order to solve the Whitham equations, it is convenient to use this form 
of the phase. 

Theorem 2.2. The phase /32, /^s) defined in \2. 6|) satisfies the following identity 
(2.26) m,/32,/33)=g(/5i,/32,/53), 

where q = q{Pi, P2, P3) is defined in i2.1!^) for X3 < and in Ii2.14\ ) for -^3 > 0, where 
Proof. Using the inverse functions f± we write p(A) and r(A) defined in ()2.5p in the form 

1 r r(x)-' ^^^-K)'^ ' 



Then the following identities hold 

r p{X)dx 

Jhui, V(A-/3i)(A-/32)(A-/33) ~ 

(2.27) ^ r f40d^ ( r + r .... i + 



24 Vife 4 V(A-/3i)(A-/32)(A-/33)(e-A) 

1 /"^^ dX f^^ ^+^^ldC 



2 7/?3 V(A-A)(A-/?2)(A-/33) A 

Using the fact that the first term in the r.h.s. of the above relation is identically zero and 
performing a change of coordinates of integration in the second term, we obtain 



1 f p{X)dX 

" hui, v/(A-/3i)(A-/32)(A-/33) 

l—f a \ I 1— At 



(2.28) 1 /' + + 



2^/27^ J-i J -I ./T^VT 



V 



2 
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For the term of containing the transmission coefficient r, we obtain for t < T 
i f T{X)d\ 



[o,,]\u/, V(A-/5i)(A-/32)(A-/33) 



2^ift V(A-A)(A-/32)(A-/33) Jps 



[2.29) = / / dfidu 



2V2n J -I j-i v^r^yr 



2 



dfidu 



2 



Using the expression ()2.28|1 and ()2.29j) and the fact that the function q{Pi, P2, Ps) defined 
in ()2.13j) is symmetric with respect to Pi, P2 and Ps, we derive the statement of the theorem 
for t < T. 

For t > T, repeating the same procedure above, we derive the relation 

. f T{\)d\ 1 /-^^ ''V-i v/A^ ^ 



(2.30) U 

^ i[o,,]\u/, V(A-A)(A-/32)(A-/33) 27r 4 ^(^ - A)(A - /32)(A - /^s) 

Using the fact that the expression ()2.28p is symmetric with respect to Pi, P2 and P^ we 
rewrite (12.281) in the form 



1 f p{\)d\ _ 1 f ^ Jf,, 



TT 



lu/. V(A-/5i)(A-/32)(A-/33) 2njp^ ^^^^ _ x){X - P2){X - P^) 



Combining the above expression with ()2.3()|1 we arrive at the expression for the phase (j) 
which coincides with q{Pi, P2, P3) defined in ()2.14|1 . □ 

In the following we show that the elliptic solution p.5|) attaches but not to the 
solution of the Hopf equation. This result is numerically obvious from Fig. El and Fig. 

Theorem 2.3. The approximate solution of the Cauchy problem / li. i)) given by formula 
h2. Sl\) in the oscillatory zone and by the Hopf solution \1.4\) outside the oscillatory zone 
is C° but not in the (x, u) plane. 

Proof. Let Uapp{x,t,e) be the r.h.s. of p.5|) and rewrite formula ()2.3ip using the Jacobi 
elliptic function dn, where 



The following identity holds |2 

1-^2 E(s) 



(ki\2zK{s)) K{s) 
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SO that, substituting the above into ()1.5p we obtain 

(2.31) Uapp{x, t,e)=/32 + /33 -/3i + 2 ~^ 

where 

(2.32) n = ^/f3i-f3^{x - 2(A +(32 + (3s)t - q). 
In the hmit (32 = (3z we have dn(z) — 1 so that 

Uapp{x,t,e) = Pi{x,t), 

where Pi{x,t) satisfies the Hopf equation because of the boundary condition ()2.9|) . 
In the limit P2 = Pi, the function dn(2;) — sechz and 

(2.33) = x- 6tps - f-m - [4t(A - Ps) - f-m + 1 r ^ySS] = 0, 

because of fll.4|) and ()2.19p . so that 

"^appiX-i 

t,e) = P3{x,t), 

where now P^ satisfies the Hopf equation because of the boundary condition ()2.10|) . There- 
fore Uappr{x, t, e) attaches C° to the solution of the Hopf equation in the limits P2 = Pi or 

P2 = P3. 

Now we show that Uappr{x,t,e) is not attached to the solution of the Hopf equation. 
For this purpose we need to evaluate dxUapp{x, t, e), namely 

f> f , ^ f,R^f,a ^ ^ . , dxPx-dA , , g^(A-/32)sn(r]/e)cn(r]/e) ^^ 

dxUappix, t, e) = dxp2 + d^Ps - d^Pi + 2—^-—— + 4 3,^^ , ^ 

dn [il/e) dn (is/^)^ 

where Q is defined in fj2.32p . For simplifying the calculation, we compute the x-derivatives 
before Ps reaches the minimum value —1. The calculation does not change in the general 
case, it is only more involved. The derivatives dxPi have been defined in ()2.23j) . We first 
consider the leading edge, namely the case P2 = Ps- Fixing P2 = v + 6, P3 = v — S, we 
obtain for 5^0, with the expansion ()2.7|1 . the relation ()2.24|1 and 

^^^2 = — ^ ,^02 ^, ^ . + 0(1), 

26{v- pi)-§^<!>{v,Pi) 

^^^3 = -— ^ ,V ^ . + 0(1). 

26{v ~ P^)-§^<l>{v,Pi) 

Furthermore as s — 0, cn{z) cos(2;), sn{z) sin(z) and dn{z) 1. Combining the 
above limits and ()2.7j) we obtain that 

dxUapp{x,t,e)\[f32=f3^] = 00, 

while dxu{x,t), where u{x,t) is the solution of the Hopf equation, remains finite at the 
leading edge. 
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At the traihng edge /?2 A or s — > 1. In this hmit cn(z) = dn{z) = sech(2;) while 
sn{z) = tanh(z) and the elliptic functions E{s) and K{s) behave as in ()2.8p . We use the 
notation P2 = v — 6,Pi=v + S and evaluate the a;-derivatives in the limit 6^0 obtaining 



OA 



6 



5 

Using the above relations, ()2.25|1 and ()2.33|1 we find 

dxUapp{x,t,e)\[p2=l3i] = +00. 

while dxu{x,t) remains finite at the trailing edge. □ 

2.2. The explicit example. In this subsection we consider the explicit example of the 
initial data 

(2.34) uoix) : ^ 



cosh (x) 

For this initial data the point of gradient catastrophe can be evaluated analytically. It is 
given by 

/S /3 

(2.35) tc=^, Xe = -^ + log((y3- 1)772), M, = -2/3. 

o / 

The increasing and decreasing part f± of the inverse function of Uq{x) take the form 



(2.36) /^(^) = lnl±^^, f'^{y) = T: 



-y 2y^TT^' 

where — 1 < y < 0. Furthermore, we obtain an analytical expression for the function 
$(A,r7) defined in (pTrjl 



(2.37) *(,,,).__L__„, 



T] — X 



^(A + 1) 



Having the above explicit expressions, we can analytically solve the system ()2.2H1 . and 
obtain the time T when (3^ = —1, 

(2.38) ^ = ^^ = -^/^' /^i = -i' /^3 = -l. 

3. Numerical solution of the KdV equation 

In this paper we are interested in the numerical solution of the KdV equation for hump- 
like initial data which are rapidly decreasing for |x| oo. Therefore it is legitimate to 
study the problem in a periodic setting of sufficiently large period. We use here a slightly 
modified version of Trefethen's code for the KdV equation (Chap. 10 in ^24J) which is 
available at 
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The basic idea of the code is the use of a discrete Fourier transform in x and an 
integrating factor such that the time derivative contains the only hnear term in u in the 
equation. Let 

u(t,k) := / u(t,x)e''''' dx, 
Jr 

be the Fourier transform of u and let be the transform of u"^. Then the transformed 
KdV equation reads 

(3.1) ut-e^ik^u + 3iku^ = 0. 
This is equivalent to 

(3.2) (e-'''''uj + 3tke-'''''u^ = 0. 

The integrating factor avoids a stiff term in the equation and thus allows for larger time 
steps. To solve equation ()3.2j) numerically we use the Fast Fourier Transform (FFT) in 
MATLAB for the x-dependence and a fourth-order Runge-Kutta method for the time 
integration. 

This code is perfectly adequate to solve the KdV equation for an e of the order 1. How- 
ever in the limit of small e we are interested in here, the aliasing error becomes significant. 
This error is due to the pollution of the numerically calculated Fourier transform u by 
higher frequencies because of the truncation of the series, see jH] for details. It becomes 
important in dealing with nonlinearities in the equations. The error can be suppressed by 
putting a certain number of the high frequency components of u equal to zero after the 
nonlinear operations. As a rule of thumb it is sufficient to put roughly 1/3 of the coeffi- 
cients equal to zero jH]. Thus effectively we are working with a lower resolution (2/3 of 
the number N of modes given below), but this avoids high frequency noise and stabilizes 
the code. 

To test the accuracy of the code we consider an exact solution to the KdV equation, the 
1-soliton u = 2/ cosh^(x — xq — 4t) for e = 1. The x-coordinate takes values in [— vr, tt]L 
where the length L is always chosen in a way that the coefficients for the initial data are 
or the order of the rounding error^ for high frequencies. This reduces the error due to the 
discontinuity of the initial data at the interval boundaries to the order of the rounding 
error. We use the 1-soliton solution at time t = with L = 10 and xq = —L as initial data 
and determine for each time step the difference of the exact and the numerical solution. 
The computation is carried out with = 2^^ modes (we will always use powers of 2 here 
since the FFT algorithm is most efficient in this case, but this is not necessary) and 4000 
time steps for t G [0,5]. The maximum of the absolute value of the difference between 
the numerical and the exact solution is shown as a function of time in Fig. 121 

An alternative test of the numerical precision is provided by conserved quantities as 
the energy, 

POO 

(3.3) E = const / {2u^ - e\l)dx. 



MATLAB works internally with a precision of 10 Due to rounding errors machine precision is 
generally limited to the order of 10^^"*. 
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t 

Figure 9. Numerical errors for the time evolution of 1-soliton initial data: 
maximum of the absolute value of the difference between exact and numer- 
ical solution (maxdiff) and deviation from energy conservation (err). 

This quantity is analytically conserved during time evolution, but numerically it will be 
a function of time due to unavoidable numerical errors. Since energy conservation is not 
implemented in the code, it provides a strong test for the numerical accuracy. We define 
the function err via 

(3.4) en -.= 1- E{t)/E{0), 

where E{t) is the numerically calculated energy which is obtained from the first component 
of the FFT of the integrand in ()3.3j) at each time step. For the example of the 1-soliton 
solution, this function is shown in Fig. IHl It can be seen that the error obtained via 
the integral quantity is typically an order of magnitude higher than the maximal local 
difference of the exact and the numerical solution. In cases where no exact solution is 
known, we will use energy conservation as an indicator of the precision of the numerical 
solution. The number of modes and the time step will be chosen in a way that the value of 
the function err is at least an order of magnitude lower than the precision of the numerical 
solution we are aiming at. 

To study the small dispersion limit of KdV solutions, we consider hump-like data of the 
shape of the 1-soliton for both signs, u = ±1/ cosh^(x). We show plots for the evolution 
of negative initial data in Fig. Eland for positive initial data in Fig. ^2 for e = 0.1 and 
t e [0, 0.4]. To show that the form of the oscillations for positive and negative initial data 
is typical, we consider the evolution of the initial data u = 2sinh(x)/ cosh^(x) at time 
t = 0.4 for e = 0.1 in Fig. [El 

For a given value of e we need a spatial resolution 27cL/N of at least e. We generally 
try to be an order of magnitude below this limit. Since we use an explicit method for 
the time integration, stability is an important issue. We find empirically that a time step 
smaller than 1/A^ leads to a stable time evolution. In Table^we give the parameters used 




o 
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Figure 10. Numerical solution of the KdV equation with e = 0.1 and 
initial data —1/ cosh^(x). 




Figure 11. Numerical solution of the KdV equation with e = 0.1 and 
initial data 1/ cosh^(x). 



in the numerical computations and the obtained numerical errors for different values of e. 
We note that the code is very efficient. Up to values of e = 10~^ the code can be run on 
standard computers without problems. 
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Figure 12. Numerical solution of the KdV equation with e = 0.1 and 
initial data 2sinh(x)/ cosh^(x) for t = 0.4. 
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logio err 
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10 


5 


4* 10-4 


-6.32 


1.25 


12 


5 


2* 10-4 


-7.79 


1.5 


12 


5 


2* 10-4 


-6.33 


1.75 


14 


5 


10-4 


-6.30 


2 


14 


5 


5* 10-*^ 


-6.29 


2.25 


16 


4 


2.5* 10-*^ 


-6.30 


2.5 


16 


4 


2.5* 10-^ 


-4.79 


2.75 


17 


4 


6.67* 10-** 


-6.16 


3 


17 


4 


6.67* 10-*^ 


-4.68 



Table 1 . Parameters used in the numerical integration of the KdV equa- 
tion for several values of e 

4. Numerical solution of the Whitham equations and of the Hope 

equation 

In this section we solve numerically the Whitham equations ()2.3|) for given initial data 
by inverting the hodograph transform ()2.11|) to obtain Pi > P2 > Ps as a function of x and 
t. Since the hodograph transform becomes degenerate at the leading and trailing edge we 
solve the system ()2.15|1 instead of ()2.1H1 to avoid convergence problems. In a similar way 
we address the implicit solution of the Hopf equation ()1.4|1 . 

Both sets of equations are of the form 

(4.1) Sii{y,},x,t) = 0, 2 = 1,...,M, 

where the Si denote some given real function of the yi and x, t. The task is to determine 
the yi in dependence of x and t. To this end we determine the yi for given x and t as the 
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zeros of the function S := Yli=i ^f- This will be done numerically by using the algorithm 
of j2ZI which is implemented as the function fminsearch in Matlab. The algorithm provides 
an iterative approach which converges in our case rapidly if the starting values are close 
enough to the solution (see below how the starting values are chosen). Generally the 
iteration is stopped for function values smaller than 10~^. For quantitative comparisons 
we calculate the zeros to the order of machine precision. 

The solution Pi{x,t) > j32{x,t) > /?3(x, t) of the Whitham equations exists for t > 
tc where tc is the time of gradient catastrophe for the solution of the Hopf equation 

Ut + 6uUx = 0. We look for a solution of the hodograph transform for t > — for a 

8 

discretized time starting with time ti close to tc- For the moment we suppose that t < T, 
where T is defined in fl2.38|) . For each fixed time, our strategy to solve ()2.11|) is the 
following: First we solve the system ()2.16p to obtain the leading edge coordinate x~{t) 
and 

At time ti the starting point for solving numerically ()2.1(jj) is given by Xq ~ Xc and 
~ Uc- Similarly we solve the equations ()2.19|1 for x~^ and Pi = (32 and /^s which fixes 
the interval x"*"]. This interval is subdivided in a number of points x„, = 1, . . . , N^- 
We use the values of the at point starting values for the iteration at point Xj. 

A typical plot is shown in Fig. |H1 Since the values of the jSt change rapidly with x near 
the leading and the trailing edge, we use a grid with one third of the points located in 
the vicinity of each of the edges and some wider grid spacing in between. To explore 
the parameter space in t and x we typically use = 30, for precision calculations we 
use = 300 and higher. Additional points for plotting purposes are obtained via cubic 
interpolation. 

For t > T, the procedure to solve the Whitham equations remains roughly the same. 
However we have to take care of the fact that the decreasing part of the initial data 
contributes to the solution of the Whitham equations when f]^ goes beyond the minimum 
value —1. Thus for t > T we determine the value xt where /^a = — 1 as a solution of 
()2.11|) . For X > we add the contributions of the decreasing part of the initial data. 
Since the algorithm j2Z| varies the values of the /3i for fixed x, one has to make sure for 
values of x near xt that the right branch of the logarithms is chosen. 

The functions (12.111) evaluated in the zero-finding algorithm contain elliptic integrals 
and integrals over the initial data which are calculated with a Chebychev collocation 
method. Elliptic integrals and functions are distributed with MATLAB where they are 
calculated with the algebro-geometric mean (see |1^) to machine precision. We use here 
the approach to hyperelliptic functions via a Chebychev collocation method presented in 
[T^ (see also Chap. 6 of [23). The elliptic integrals and functions can be calculated to 
machine precision with this approach which is checked via internal tests and a comparison 
with the functions calculated with MATLAB. The reasons for the use of the Chebychev 
method are twofold. First we can apply similar techniques to calculate the non-standard 
integrals in 1)2.111) with machine precision. And secondly we develop in this way an 
approach to the one-phase Whitham equation which is directly open to a generalization 
to the multi-phase Whitham equations since the Chebychev code can handle hyperelliptic 
Riemann surfaces of in principle arbitrary genus. 
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Let us briefly summarize the Chebychev approach, for details see [13 1211 El CSI • The 
Chebyshev polynomials T„(x) are defined on the interval I = [—1, 1] by the relation 

Tnicosit)) = cos(nt) , where x = cos(t) , t G [0,7r] . 

The addition theorems for sine and cosine imply the recursion relation 

(4.2) - = 2T„(x) 

n + 1 n — 1 

for their derivatives. The Chebyshev polynomials are orthogonal on / with respect to the 
hermitian inner product 

dx 

U.9)= / f{x)-g{x) 



2 



1 — X 

We have 

vr 

(4.3) (TjYi, 7n) 2" 

where Cq = 2 and Q = 1 otherwise. A function / on / is approximated via Chebychev 
polynomials, / ~ S^=o '^n^n(a;) where the spectral coefficients a„ are obtained by the 
conditions f{xi) = ^n=o^"-'^-n{xi), I = 0,...,N. This approach is called a collocation 
method. If the collocation points are chosen to be xi = cos{ttI/N), the spectral coefficients 
follow from / via a Discrete Cosine Transform for which fast algorithms exist. 

The fact that / is approximated globally by a finite sum of polynomials allows us to 
express any operation applied to / approximately in terms of the coefficients. Let us 
illustrate this in the case of integration. So we assume that f = Pn = X]n=o '^nTn and we 
want to find an approximation of the integral for p^, i.e., the function 



F{x) = j J{s)ds 



so that F'{x) = f{x). We make the ansatz F{x) = J2n=o ^nTn{x) and obtain the equation 

N N 

F' = Y,bnZ = J2<^nTn = f. 

n=0 n=0 

Expressing T„ in terms of the via ()4.2|1 and comparing coefficients we get the equations 

h = , On = ^ for < n < A/ , = —r- . 

between the coefficients which determine all bi in terms of the a„ except for bo- This 
free constant is determined by the requirement that F{—1) = which implies (because 
T„,(-l) = (-1)") 

N 

bo = -J2^-irbn . 

n=l 

From the coefficients bn we can also find an approximation to the definite integral J^-^ f{s) ds 
F(l) by evaluating 

N [N/2\ 
gAr(l) = ^bn = 2 ^ 62^+1 • 
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Thus to find an approximation of tlie integral of a function / we proceed as described 
above, first computing tlie coefficients a„ of /, computing tlie 6„ and tlien calculating the 
sum of the odd coefficients. 

There are two types of integrals in ()2.1H) . Integrals of the form 



are calculated with the Chebychev integration routine after the transformation n = a + 
{b — a)(l + uY where y & I. The reason for this transformation is to obtain a smooth 
integrand free of square roots which is important for an efficient use of spectral methods. 
After the square root substitution we map the integration to the interval / with a linear 
transformation. The second type of integral is of the form 



This integral is calculated by expanding the function / in terms of Chebychev polynomials 
and by using the orthogonality relation ()4.Hj) . The precision of the numerical calculation is 
controlled via test integrals of known functions and via the spectral coefficients: the latter 
have to be of the order of machine precision for n ~ to provide sufficient resolution. 
In our examples 32 to 128 polynomials were always sufficient to fulfill this requirement. 
Note that there are problems with the integrands if (3^ ~ —1 since the f± diverge there. 
This would require an additional coordinate transformation to obtain a smooth integrand 
necessary for an efficient spectral approximation. Here we are, however, able to perform 
the required integration by hand. This has the additional benefit to provide a faster code, 
but the method is able to handle general initial data of the considered form. 

To obtain a solution of the Hopf equation we choose again a convenient numerical grid 
in t and x G [xl,xr\, where xl and xr are x- values where the solution is single valued. 
We solve equation ()1.4|) for fixed t and all corresponding values of x starting from xl to 
x~ with starting value ^ = x as described above. Successively we repeat the procedure 
for greater values of Xi on the grid with starting value the value of ^ found in the 
preceding step. We find that this approach works well even for x = x~ very close to the 
point of gradient catastrophe Xc- Similarly we solve the Hopf equation starting from xr 
to x"*", the x-coordinate of the trailing edge. The solution of the Hopf equation in the 
region where it is multivalued, is obtained in Fig. Elby plotting the contour of zero values 
of the function x — Qtu ± ln((l + ^Jl + u)/ y/—u). It can be seen form Fig. ^Jthat the 
matching of the solution to the Whitham equations to the solution of the Hopf equation 
is always C^, but the solution of the latter in the multivalued region does not coincide 
with the former. The Whitham zone grows in time as can be seen from Fig. [7| 

To sum up we have shown in this section that we are able to obtain the numerical solu- 
tion of the Whitham equations and the asymptotic approximation to the small dispersion 
limit of the KdV equation with machine precision. The given method is in principle open 
to deal with general hump-like data with a single minimum. A generalization to a higher 
genus situation appears to be straight forward. Numerical solutions of the Whitham 
equations in the genus two case has been obtained in |22I for the Benjamin-One equation. 



(4.4) 




(4.5) 
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Figure 13. Solution of the Hopf equation and the Whitham equation 
(thick hne) for the initial data uo{x) = — 1/ cosh^x at various times. 

5. Comparison of the small dispersion KdV solution and its 

approximation 

In the previous sections we have shown how to obtain numerical solutions to the KdV 
equation in the small dispersion limit as well as the asymptotic solution ()1.5|1 which 
follows from the solution of the Whitham equations, both with controlled precision. This 
enables us to present a quantitative comparison of the KdV solution and the asymptotic 
approximation for several values of the dispersion parameter e. The code works for general 
rapidly decreasing initial data with a single hump, for positive initial data, see Fig. 
This shows numerically that the formula for the phase in Theorem 2.2, which was originally 
obtained only for negative initial data, can also be used for positive initial data. 

In the following we will study as an example the initial data 

Uo{x) = 

cosh X 

The quality of the numerics allows to reach values of e of the order of 10^^ without 
problems. For small e the asymptotic approximation is almost identical to the KdV 
solution in the Whitham zone as can be seen in Fig. |31 Thus for small e it makes little 
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(a) KdV solution (b) Asymptotic solution 

Figure 14. In (a) the numerical solution of KdV for the positive initial 
data Uo{x) = 1/ cosh^x and in (b) the corresponding asymptotic formula 
(fT3|l and jrij) are plotted for t = 0.35 and e = 10"^ 



sense to plot both solutions in one figure as in Fig. |2lfor e = 10~^. We show the envelope 
of the asymptotic solution and the solution ()1.4|1 together with the solution of the 
KdV equation in Fig. ^1 The shape of the envelope follows from ()2.3H1 to be given by 
P1-P2 + P3 and P1 + P2- Ps- 

The numerical precision enables us to study quantitatively the difference between the 
KdV solution and the asymptotic solutions p.5|) and ()1.4|) and to establish where the 
approximation is satisfactory and where it is not. We show this difference for various 
values of e in Fig. ^Jat time t = 0.4. 

From Fig. ^Jand Fig.lHl it is clear that the error is decreasing with e. Indeed in the first 
plot, the highest peak is of the order of 0.15, while in the last plot it is of the order 0.06. 
It is also obvious that the asymptotic formula p.5|l gives a satisfactory description of the 
oscillations within an interval in the Whitham region [x~{t),x~^{t)]. At the boundaries of 
this zone, the highest peaks in the difference appear. 

A similar behavior can be observed for the asymptotic solution ()1.4|) . The smaller e, the 
better the approximation, which is always worst at the boundary of the Whitham zone. 
It is also clear that there are in any case oscillations of the KdV solution for values of 
X < x~ (t) for t > tc whereas the solution p.4|) of the Hopf equation, which is considered 
as an approximation there, obviously does not show any oscillations. As can be seen from 
Fig. El this oscillatory zone shrinks when e becomes smaller. For values of x > x~^{t), 
no oscillations are observed but the difference between the KdV solution and the solution 
()1.4|1 is biggest at the boundary of the Whitham zone and goes asymptotically to zero. 
The zone where the solutions differ considerably also shrinks with decreasing e. Below we 
will study these features in more detail. 
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Figure 15. Solution of the KdV equation with initial data uq = 
— 1/ cosh(a;)^ for t = 0.4 together with the solution ()1.4j) and the envelope 
of the solution jUni) for e = 0.01. 



We will present certain characteristic quantities as the difference between the exact and 
the approximative solution in dependence of e and show loglog-plots of these quantities. 
If appropriate we perform a linear regression analysis to identify a scaling behavior in e 
of the studied quantity. We briefly summarize the basic relations for a linear regression 
(see e.g. Given a set of real points i/i, Zi with i = 1, .., M, we perform a least square 
fitting to the line y = az + b where with {z, y are the mean values) 

M M M 

(5.1) s^^ = ^{zi - zf, s^y = ^{zi - z){yi - y), Syy = ^{yi - 2/)^ 

i=l i=l i=l 

the parameters a and b are given by 

(5.2) b = ^, a = y-bz. 



s 



zz 



For the correlation coefficient being a measure of the quality of the fitting, and the for 
the standard errors of a and b one has 




We will only present the results of linear regression if the correlation coefficient is at 
least of the order of 0.99. Even in these cases the results have to be taken with care. 
We only use a small number of points, but more importantly, we only study values of 
e between 10~^ and 10^^. Thus these plots show a scaling law in this regime. Further 
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Figure 16. The blue line describes the difference between the numerical 
solution of the KdV equation and the asymptotic formula (jl.Sj) for the 
initial data uq{x) = — 1/ cosh^x and for t = 0.4. The green lines represent 
the difference between the numerical solution of the KdV equation and the 
Hopf solution 



analytical work will have to show whether these results hold for more general values of e, 
and what are the precise values of the coefficients. 

Whitham zone. From Fig. ITT)1 and Fig. IHlit is obvious that the asymptotic approximation 
does not give a satisfactory description of the KdV solution close to the boundary. 
However one can identify an interior zone where this approximation gives a rather accurate 
description of the KdV zone. There is an obvious arbitrariness in the definition of such 
an interior zone since the difference between KdV solution and approximation is not 
constant there. As can be seen from Fig. [T7I there is also an error in the phase. Close to 
X = —2.5 the KdV solution and the asymptotic solution are in phase, and consequently 
the difference is minimal there. 

A possible definition of the interior zone is simply to cut off the big oscillations at the 
edges. If one cuts off just one full oscillation, the cut off zone scales roughly like e in 
accordance with the fact that there are oscillations of the order 1/e in the Whitham zone. 
For the case t = 0.4 and e = 10~^, we obtain Fig. El 



SMALL DISPERSION LIMIT OF KDV 



29 




X 



Figure 17. Difference between the numerical solution of the KdV equation 
and the asymptotic formulas ()1.5|) for t = 0.4 and e = 10~^ in the 'interior' 
Whitham zone. 



To define an error (err^id) in the interior we take the maximum of the absolute value of 
the difference between the solutions in the vicinity of the center of the zone {x ~ —2.66). 
Notice that there is a certain crudeness in this definition of the error since this maximum 
will occur for different e at slightly different values of x. Because of the error in the phase 
there are large differences in the error for different e if one takes the difference at the same 
X value in all cases. The so defined error is shown in Fig. ^1 where it can be seen that it 
decreases almost linearly with e. Linear regression analysis yields a = 1.0049, b = 0.1005 
with (Ta = 0.05, (Tfe = 0.024, and a correlation coefficient r = 0.998. 

The maximal error near the boundaries of the Whitham zone is shown in Fig. ^1 There 
seems to be no obvious scaling of these errors with e which is probably due to the rapid 
oscillations of the error because of the error in the phase. 

Oscillatory zone with x < x~{t). As can be seen in Fig. Eland Fig. there are always 
oscillations of the KdV solution for t > tc which do not occur for the asymptotic solution 
fll.4|l . We show this region in detail in Fig. 1201 

It can be seen that the oscillatory zone clearly shrinks with decreasing e. There is no 
obvious rigorous definition of the end of the oscillatory zone. The difference will eventually 
be of the order of the numerical error for the KdV solution. We define the end of the 
oscillatory zone as the x-value where the amplitude of the difference of the solutions is 
smaller than 10~^ which is of the order of the numerical accuracy we have used. The 
width of this zone, ^^opf •= ^hopf/^~ ~ 1 in dependence of e is shown in Fig. EH It can be 
seen that the oscillatory zone shrinks to zero with e. Asymptotically the boundary of this 
zone approaches the boundary of the Whitham zone. We find that Ay^^pf scales roughly 
as e^/^. The linear regression analysis yields a = 0.761 and b = —0.789 with cTq = 0.028, 
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Figure 18. Maximum of the difference between the KdV numerical so- 
lution and the asymptotic formula ()1.5|1 in the interior Whitham zone for 
t = 0.4 in dependence of e; the data can be fitted by a straight line y = az+b 
with a = 1.0049 and b = 0.1005. 
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Figure 19. Maximum of the difference between the KdV numerical solu- 
tion and the asymptotic formula (jl.Sp at the boundaries of Whitham zone 
for t = 0.4 in dependence of e (* near x~ and o near x^). 




ab = 0.013 and r = 0.999. However, this result has to be taken with care in view of the 
low number of points and the arbitrariness in the definition of the zone. The error err^^^^ 
in this zone is measured by the maximum of the absolute value of the difference of the 
amplitude between the KdV solution and the Hopf solution. This maximum value always 
occurs close to the boundary of the Whitham zone. The error err^^^^ is shown in Fig. |221 
It decreases roughly as e^/^. Linear regression analysis yields a = 0.346, b = 0.250 with 
aa = 0.025, ab = 0.012 and r = 0.996. 

Zone with x > x^{t). A similar behavior is found in the zone for x > x~^{t) with the 
exception that there are no oscillations. The difference of the asymptotic solution ()1.4|) 
and the KdV solution is biggest for x~^{t) and decreases monotonically to the order of the 
numerical precision. 

We define the width of the zone as the region where the absolute value of the difference 
is bigger than 10""^. The values of the quantity ^topf ■~ ^~^hopf/^~^ shown in Fig. 1^ 
in dependence on e. It can be seen that this zone shrinks with e to zero, i.e., its boundary 
approaches asymptotically the boundary of the Whitham zone. Apparently there is no 
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Figure 21. The size of the left oscillatory zone ^hopf ^ function of e 
for t = 0.4 in dependence of e; the data can be fitted by a straight line 
y = az + b with a = 0.761 and b = —0.789. 
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Figure 22. Maximum of the absolute value of the difference err^^^^ be- 
tween the numerical solution of the KdV equation and the Hopf solution 
()1.4|1 for t = 0.4 in dependence of e; the data can be fitted by a straight line 
y = az + b with a = 0.346 and b = 0.25. 
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Figure 23. Difference between the numerical solution of the KdV equation 
and the Hopf solution ()1.4|1 for t = 0.4 and x > x^(t). 



simple e-dependence of this quantity, at least not for the number of points used in the 
plot. 

The maximal error in this zone always occurs at the boundary to the Whitham zone. 
As can be seen in Fig. |2S1 this error decreases roughly as ^/e. Linear regression analysis 
yields a = 0.525, b = 0.554 with = 0.017, at = 0.008 and r = 0.999. 
Breakup time. As discussed above the asymptotic approximation of the small dispersion 
KdV solution is always worst near the boundary of the Whitham zone. At the breakup 
time tc being defined as the time of gradient catastrophe of the solution to the Hopf 
equation, the Whitham zone consists only of one point. In Fig. El it can be seen that the 
solution to the KdV equation always forms oscillation for times t < tc- Thus the first 
oscillation in the Whitham zone gives only a very crude approximation to the oscillation 
of the KdV solution in this zone. In addition there are typically several oscillation of 
the KdV solution outside the Whitham zone at this time. The discrepancy between the 
asymptotic and the KdV solution is quite pronounced even for small e as can be seen from 
Fig. EElfor e = 10~^. The KdV solution develops roughly the same number of oscillations 
before tc as in Fig. IHl but on a smaller scale. As already stated, the oscillatory zone 
shrinks with decreasing e, but close to the breakup the difference between the asymptotic 
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Figure 24. The size of the zone A^^pj els a function of e for t = 0.4. 




Figure 25. The maximal error in the zone x > x'*' as a function of e for 
t = 0.4; the data can be fitted by a straight hne y = az + b with a = 0.525 
and b = 0.554. 



solutions and the KdV solution remains considerable, even for small e. As can be seen in 
Fig. 123 the difference between asymptotic solution and KdV solution is biggest close to 
the breakup time. 
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t=0.216 t=0.218 




Figure 26. KdV solution and asymptotic solution for e = 10 ^ close to 
the breakup time. 

Time dependence. The qualitative behavior of the difference of the small dispersion KdV 
solution and the asymptotic solution as outlined above is characteristic for all times as can 
be seen from Fig. |2ZI The difference is always biggest at the boundary of the Whitham 
zone. The absolute value of this difference is almost constant in time at the boundary. 
Therefore the discrepancy between the approximation and the solution is biggest close 
to breakup. The interior zone where the asymptotic solution gives a very satisfactory 
approximation grows with time. 

6. Outlook 

In this paper we have presented a quantitative numerical treatment of the solution to 
the KdV equation in the small dispersion hmit for initial data with compact support and a 
single hump. The same has been achieved for the asymptotic formulas for the same initial 
data. The code for the KdV solution is set up to handle general initial data with compact 
support, thus the inclusion of initial data which develop multi-phase solutions is directly 
possible. The approach to the solution of the one-phase Whitham equations is also open 
to a generalization to multi-phase case. In particular the existence and uniqueness of the 
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Figure 27. Difference of the KdV solution and tlie asymptotic solution 
for e = 10~^ for different times. 



solution of the two-phase Whitham equations has been obtained in |17j for generic initial 
data where /_ has negative fifth derivative. Since the Chebychev code [1211211 to calculate 
the characteristic quantities on the elliptic surface in the present paper is able to deal with 
almost arbitrary hyperelliptic surfaces, the corresponding models can be studied. This 
will be the subject of future work. 

The results of the previous section show that the asymptotic approximation is worst at 
break up point and at the boundary of the Whitham zone. According to the conjecture 
in jin],|^ the solution of the KdV equation, near the point of gradient catastrophe, is 
well described by 

{6.1) u{x,t,e) = u,+ i^-j ^(--(^j {x - X, - 6u,{t - Q); - i^-j {t - Q j , 
where 
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and the function F{X; T) is the solution of the fourth order ODE of the Painleve I 
hierarchy 

(6.2) -QTF + F^ + FF" + 1 /2(F')^ + 1/lOF"" = X. 

The above equation first appeared in the double scaling limits of one matrix models for 
the multicritical index m = 3 and for T = P|. It is conjectured that the equation ()fi.2j) 
has a smooth real solution on the real line with asymptotic behavior F{X, 0) = ztXa for 
X ±oo. We will investigate numerically the asymptotic description given by ()6.H) in a 
subsequent publication. 

The description of the left oscillatory zone outside the Whitham zone, should be ob- 
tained in the spirit of matrix models, performing a sort of double scaling limit. According 
to the ansatz in the envelope of the oscillations is given by the Hastings McLeod 
solution of the second Painleve equation. This result is in accordance with the rigorous 
results obtained in the double scaling limit in one-matrix models |ni)0;IZl- 

Regarding the right boundary of the Whitham zone, at the moment, there is no theo- 
retical ansatz for an asymptotic description. Our numerical results could possibly provide 
a hint for an ansatz for the asymptotic behavior. 

It would also be very interesting to perform similar numerical investigations for the 
semiclassical limits of the focusing |2S1 and defocusing |2I] nonlinear Schrodinger equation. 
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